#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*
 * 幂法求矩阵的模最大特征值和对应特征向量
 * 参数一：矩阵
 * 参数二：矩阵的行数
 * 参数三：矩阵的列数
 * 参数四：所要求的精度
 * 参数五：传出参数，传出矩阵的特征值
 * 返回值：矩阵特征值对应特征向量所在数组的首地址
 */
double *PowerLaw(double Matrix[3][3],int line,int column,double precision,double *Eigenvalue);
/*
 * 幂法的子函数，用于找数组中最大的值并返回
 */
double FindMax(const double *Array,int num);

int main() {

    double Matrix01[3][3] = {
            {2,3,2},
            {10,3,4},
            {3,6,1}
    };

    double Matrix02[3][3] = {
            {2,4,6},
            {3,9,15},
            {4,16,36}
    };

    double Matrix03[3][3] = {
            {1,-3,2},
            {4,4,-1},
            {6,3,5}
    };

    double Matrix04[3][3] = {
            {0.7312,1,0},
            {1,1.7321,0},
            {0,1,2.7321}
    };

    double Matrix05[3][3] = {
            {2,1,0},
            {1,3,1},
            {0,1,4}
    };
    double Eigenvalue;
    double *Ans = PowerLaw(Matrix05,3,3,0.00001,&Eigenvalue);



    printf("Hello, World!\n");
    return 0;
}

double *PowerLaw(double Matrix[3][3],int line,int column,double precision,double *Eigenvalue){
    double PreAns[line];    //用于保存上一次计算时的答案，最终用于计算精度
    double Yk[line];    //用于保存每次Xk与矩阵计算后的结果
    double *NowAns = malloc(sizeof(double )*line);  //用于保存每次的Xk向量的值
    memset(PreAns,0, sizeof(double )*line); //设置preAns的初值

    //用于设置Xk初向量
    for (int i = 0; i < line; ++i) {
        NowAns[i]=1;
    }

    int flag=1;     //精度标记
    while (flag){

        //计算矩阵乘以初向量
        memset(Yk,0, sizeof(double )*line);
        for (int i = 0; i < line; ++i) {
            for (int j = 0; j < line; ++j)
                Yk[i]+=Matrix[i][j]*NowAns[j];  //矩阵与列向量乘法
        }

        *Eigenvalue = FindMax(Yk,line); //获取当前Mk值

        for (int i = 0; i < line; ++i){     //更新当前Xk向量值
            PreAns[i]=NowAns[i];    //保存当前xk值，用于精度计算
            NowAns[i]=Yk[i]/(*Eigenvalue);  //将Yk向量除以其最大分量的值
        }

        //查看当前Yk值Mk值和Xk值
        for (int i = 0; i < line; ++i) {
            printf("Yk%d: %.4lf\t",i+1,Yk[i]);
        }
        for (int i = 0; i < line; ++i) {
            if(i==0)
                printf("E: %.4lf\t",*Eigenvalue);
            printf("Xk%d: %.4lf\t",i+1,NowAns[i]);
        }
        printf("\n");

        //确定当前Xk的精度是否符合要求
        for (int i = 0; i < line; ++i) {
            if(fabs(NowAns[i]-PreAns[i])>precision)
                break;
            if(i == line-1)
                flag=0;
        }
    }

    return NowAns;
}

double FindMax(const double *Array,int num){
    double Max = *(Array);
    for (int i = 0; i < num; ++i) {
        if(*(Array+i)>Max)
            Max=*(Array+i);
    }
    return Max;
}